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Preface 


The  work  reported  herein  was  conducted  as  part  of  the  Water  Quality 
Research  Program  (WQRP),  Work  Unit  32808,  “Model  of  Contaminant 
Transport  and  Fate  at  Corps  Projects.”  The  WQRP  is  sponsored  by  the 
Headquarters,  U.S.  Army  Corps  of  Engineers  (HQUSACE),  and  is  assigned  to 
the  U.S.  Army  Engineer  Waterways  Experiment  Station  (WES)  under  the 
purview  of  the  Environmental  Laboratory  (EL).  Funding  was  provided  under 
Department  of  the  Army  Appropriation  96X3121,  General  Investigation.  The 
WQRP  is  managed  under  the  Environmental  Modeling,  Simulation,  and 
Assessment  Center  (EMSAC),  Dr.  John  W.  Barko,  Director  for  EL. 

Mr.  Robert  C.  Gunkel  was  Assistant  Manager  for  the  WQRP.  Program  Monitor 
during  this  study  was  Mr.  Frederick  B.  Juhle,  HQUSACE. 

The  Principal  Investigator  of  Work  Unit  32808  was  Dr.  Mark  S.  Dortch, 
Chief,  Water  Quality  and  Contaminant  Modeling  Branch  (WQCMB), 
Environmental  Processes  and  Effects  Division  (EPED),  EL.  The  work  reported 
herein  was  conducted  by  Dr.  Raymond  S.  Chapman,  Ray  Chapman  and 
Associates,  Vicksburg,  MS,  and  Mr.  Terry  K.  Gerald,  AScI,  Inc.,  Vicksburg, 

MS,  under  contract  to  WES.  This  report  was  prepared  by  Dr.  Chapman  and 
Mr.  Gerald.  Dr.  Dortch  monitored  the  contract  and  reviewed,  edited,  and  revised 
the  report. 

This  study  was  conducted  under  the  general  supervision  of  Dr.  Richard  E. 
Price,  Acting  Chief,  EPED,  and  Dr.  John  Harrison,  Director,  EL.  This  report 
was  reviewed  by  Messrs.  Thomas  Cole  and  Ross  Hall,  WQCMB. 

At  the  time  of  publication  of  this  report.  Director  of  WES  was  Dr.  Robert  W. 
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ICM,  ’  Miscellaneous  Paper  W-97-1,  U.S.  Army  Engineer  Waterways 
Experiment  Station,  Vicksburg,  MS. 


The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication,  or 
promotional  purposes.  Citation  of  trade  names  does  not  constitute  an  official 
endorsement  or  approval  of  the  use  of  such  commercial  products. 


1  Introduction 


Background 

CE-QUAL-ICM  (Cerco  and  Cole  1995)  is  a  three-dimensional  (3-D)  water 
quality  model  based  on  the  finite  volume  approach.  There  are  two  versions  of 
this  model,  a  eutrophication  version  documented  by  Cerco  and  Cole  (1995),  and 
a  toxic  substance  version  referred  to  as  ICMATOXI  documented  by  Wang  et  al. 
(1996).  The  project  that  funded  the  development  of  ICM/TOXI  also  funded  the 
present  study  reported  herein. 

Neither  version  of  CE-QUAL-ICM  (ICM)  solves  for  hydrodynamics,  so  this 
information  must  be  supplied  by  a  hydrodynamic  model  for  driving  transport  in 
the  ICM  model.  The  ICM  model  is  usually  linked  to  output  firom  CH3D-WES 
(Johnson  et  al.  1991),  a  3-D  finite  difference  hydrodynamic  model  based  on  a 
structured  grid  scheme.  Structured  grids  assume  that  computational  cells  are 
ordered  along  columns  and  rows,  whereas  unstructured  grids  allow  cells  to  be 
arranged  in  a  nonordered  fashion.  Thus,  unstructured  grids  allow  greater 
flexibility  for  describing  the  geometry  of  the  domain. 

RMAIO  (Norton,  King,  and  Orlob  1973;  Thomas  and  McAnally  1990)  is  a 
3-D  finite  element  hydrodynamic  model  frequently  used  by  the  U.S.  Army 
Engineer  Waterways  Experiment  Station  (WES).  RMAIO  is  based  on  an 
unstructured  grid  approach.  There  was  a  need  for  developing  a  linkage  to  the 
RMAIO  model  to  allow  greater  flexibility  provided  by  unstructured  grids  for 
conducting  water  quality  and  contaminant  model  studies. 


Objective  and  Scope 

The  objective  of  this  study  was  to  develop  software  to  provide  linkage  of 
RMAIO  output  to  the  ICM  code  and  to  test  the  success  of  the  linkage  for  single 
test  cases.  The  software  development  consisted  of  basically  three  parts: 

(a)  development  of  a  code  (MAPPER)  to  map  the  finite  element  grid 
configuration  and  geometry  information  into  a  file  that  can  be  interpreted  by  the 
ICM  code;  (b)  development  of  an  RMAIO  postprocessor  code  (FEMCONVT)  to 
convert  RMAIO  velocity  and  water  surface  information  at  nodes  to  flows  and 
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cell  areas  and  volumes  for  each  ICM  cell;  and  (c)  modifications  within  ICM  to 
conform  to  the  RMAIO  linkage.  Linkage  testing  consisted  of  checking  local  and 
global  volume  and  mass  conservation  and  examining  the  transport  against  known 
solutions.  It  is  important  to  maintain  mass  conservation  in  water  quality  models 
since  they  are  based  on  the  mass  conservation  principle. 

This  report  is  organized  into  chapters  on  Introduction,  Implementation, 
Testing,  and  Summary  and  Recommendations.  Appendixes  A-C  deal  with  the 
MAPPER,  FEMCONVT,  and  ICM  codes. 


2 
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2  Implementation 


The  basic  requirements  for  the  design  and  development  of  an  unstructured 
grid  linkage  methodology  for  ICM  are  as  follows;  (a)  generating  a  geometry  file 
that  relates  the  elements  and  their  geometric  attributes  for  an  unstructured  grid 
hydrodynamic  model,  such  as  RMAIO,  to  flow  faces  and  computational  cells  of 
ICM;  (b)  assigning  a  tagging  system  to  relate  flows  from  the  hydrodynamic 
model  to  ICM  flow  faces;  (c)  establishing  a  convention  for  positive/negative 
flow  directions;  and  (d)  computing  integrated  flows  across  and  normal  to  each 
hydrodynamic  element  face.  The  linkage  program  MAPPER  conducts  the  first 
three  tasks,  while  the  linkage  program  FEMCONVT  conducts  the  last. 


MAPPER 

The  linkage  geometry  and  map  files  are  derived  from  the  RMAIO  element 
connection  file  that  is  generated  for  a  hydrodynamic  run.  A  C  language  program 
(MAPPER.C)  was  written  to  generate  this  information.  MAPPER.C  is  described 
and  listed  in  Appendix  A.  Descriptions  of  the  RMAIO  node  and  element 
connectivity  scheme  were  obtained  from  the  RMAIO  documentation  (Thomas 
and  McAnally  1990;  Brigham  Young  University  1994).  The  element  connection 
file,  generated  by  an  RMAIO  preprocessor  code,  which  resides  within  the  WES 
Coastal  and  Hydraulics  Laboratory  (CHL),  provides  element  and  node 
numbering,  the  3-D  coordinates  of  each  node  within  the  grid,  and  boundary 
designation  information.  This  file  is  read  by  MAPPER.C,  which  generates  the 
geometry  (GEO)  and  map  (MAP)  files  required  as  input  to  ICM. 

To  accommodate  unstructured  grid  hydrodynamics,  the  GEO  file  was 
modified  so  that  a  one-dimensional  array  of  centroid-based  grid  distances,  which 
are  direction  sensitive,  is  computed  and  replaces  the  original  box  lengths  ( i.e., 
BL(IB(F),QD(F))).  The  centroid-based  grid  lengths  data  consist  of  (a)  distance 
from  centroid  of  IB  cell  to  flow  face  (BID);  (b)  centroid  to  centroid  distance 
between  IB  and  JB  cells  (BI);  (c)  centroid  to  centroid  distance  between  the  ILB 
cell  and  IB  (BIL);  and  (d)  centroid  to  centroid  distance  between  the  JB  cell  and 
JRB  cell  (BIR)  (see  Figure  1).  The  specification  of  the  ILB,  IB,  JB,  and  JRB 
ICM  cell  numbers  (Figure  1)  is  accomplished  through  examination  of  the 
element  connection  file.  Specifically,  when  elements  that  share  common  mid¬ 
side  nodes  are  identified,  a  cell  face  between  boxes  IB  and  JB  is  defined.  Cells 
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FF 


Figure  1 .  Grid  nomenclature 

to  the  right  or  left  of  cells  with  a  common  flow  face  are  determined  by  selecting 
the  cell  with  the  shortest  centroid  distance.  For  example,  as  shown  in  Figure  1, 
the  centroid  distance  firom  Cells  JB  to  JRB  is  shorter  than  the  distance  from 
Cells  JB  to  Jl;  thus.  Cell  JRB  is  selected  as  the  adjacent  right  cell,  or  the  far 
downstream  cell  for  the  QUICKEST  weighting,  and  BIR  is  the  distance  from  JB 
to  JRB. 

Due  to  the  fact  that  unstmctured  grids  have  no  preferred  grid  direction,  a 
convention  has  been  adopted.  A  flow  direction  for  the  establishment  of  the  GEO 
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and  MAP  files  is  assumed  based  upon  the  cell  numbering  scheme  obtained  from 
RMAIO.  In  the  horizontal  direction,  flow  is  defined  to  be  positive  when  it  passes 
from  a  smaller  RMAIO  element  number  to  a  larger  RMAIO  element  number.  In 
the  vertical  direction,  flows  are  defined  positive  upward,  which  is  the  convention 
used  in  both  RMAIO  and  ICM.  When  the  flow  field  is  read  in  from  the  HYD 
file,  the  actual  flow  directions  relative  to  the  assumed  flow  directions  are 
determined  so  that  the  proper  cells  are  sampled  for  the  QUICKEST  advection 
operators. 

The  general  concept  of  the  MAP  file  was  unaltered;  however,  additional 
information  was  required  to  relate  the  RMAIO  and  ICM  grids.  The  processor 
program  MAPPER  assigns  ICM  cell  numbers  to  RMAIO  element  numbers  and 
writes  the  ordered  pairs  as  the  first  set  of  records  in  the  MAP  file.  The 
correspondence  between  hydrodynamic  elements  and  ICM  cells  is  arbitrary  in 
that  the  connection  is  accomphshed  via  flow  face  definitions.  This  approach  is 
adopted  so  that  the  existing  convention  for  numbering  cells  and  flow  faces  in 
ICM  is  unaltered. 

The  RMAIO  element  edges  and  ICM  flow  faces  are  related  in  the  MAP  file 
via  a  tag.  In  the  horizontal  direction,  the  flow  faces  between  adjacent  cells  in 
ICM  are  tagged  to  mid-side  nodes  in  RMAIO,  which  are  defined  along  element 
edges.  Vertical  flow  faces  in  ICM  are  tagged  to  the  top  of  each  element  in 
RMAIO.  These  tags  are  used  via  translation  arrays  to  relate  RMAIO  flow 
locations  to  the  corresponding  sequentially  numbered  ICM  flow  face  numbers. 

Boundaries  within  the  RMAIO  grid  are  specified  via  nodal  codes  that  define 
land,  inflow,  and  tidal  or  head  boundaries.  These  codes  are  used  in  MAPPER  to 
specify  boundaries  in  the  MAP  file. 


FEMCONVT 

The  Fortran  language  program  FEMCONVT.f  was  written  by  Dr.  R.  C. 
Berger  of  the  WES  CHL  to  read  a  binary  RMAIO  hydrodynamic  output  file  and 
convert  nodal  velocities  and  surface  elevations  into  cell  face  areas  (including 
planar  surface  areas),  volumes,  and  flows,  which  are  written  to  another  output 
file  (HYD  file).  A  listing  of  FEMCONVT  is  provided  in  Appendix  B.  The 
output  of  FEMCONVT  is  segregated  into  horizontal  and  vertical  con^Mjnents.  In 
the  horizontal  direction,  the  information  provided  is  time-varying  horizontal  flow 
face  areas  and  flows  for  each  mid-side  node.  The  vertical  part  of  the 
hydrodynamic  output  provides  element  number,  element  planar  surface  area, 
volume,  and  the  vertical  flow  defined  on  the  top  face  of  each  element. 


ICM  Modifications 

The  original  advection  operator  in  ICM  was  designed  to  work  with  structured 
grid  hydrod5mamic  models  such  as  CH3D-WES  (Chapman  1988).  As  a  result, 
the  grid  is  described  by  rows  and  columns  of  grid  cells  and  box  lengths.  These 
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box  lengths  are  used  to  compute  the  QUICKEST  multipliers  used  for  transport. 
Due  to  the  success  realized  during  numerous  applications  of  the  QUICKEST 
transport  algorithm  in  both  the  ICM  and  CE-QUAL-W2  models  (Dortch, 
Chapman,  and  Abt  1991;  Chapman  1992;  Chapman  and  Cole  1992),  a  similar 
approach  was  adopted  for  the  unstructured  version  of  ICM.  The  unstructured 
grid  modifications  made  to  the  original  version  of  ICM  are  presented  in 
Appendix  C.  For  ICM  to  accommodate  imstructured  hydrodynamic  models  such 
as  RMAIO,  the  QUICKEST  multipliers  had  to  be  recast  in  terms  of  element 
centroid  to  centroid  distances.  As  previously  discussed,  the  linkage  processor 
program  MAPPER  generates  the  four  length  scales  needed  to  develop  the 
QUICKEST  multipliers.  Using  these  length  scales,  the  QUICKEST  multipliers 
were  rewritten  and  checked  to  ensure  that  the  advection  multipliers  summed  to 
one,  and  the  diffusive  multipliers  summed  to  zero. 


Linkage  Procedure 

To  complete  a  linkage  of  RMAIO  and  ICM,  the  steps  outlined  below  must  be 
completed. 

a.  A  binary  output  file  generated  by  RMAIO,  which  contains  the  time- 
invariant  grid  and  time- varying  hydrodynamic  data,  must  be  obtained 
from  the  hydrodynamic  modeling  team.  The  name  of  this  file  is  specified 
as  RMAIO  input,  so  it  is  arbitrary. 

b.  With  regard  to  the  time-invariant  grid  data,  a  program  named  R4ICR10.f, 
written  and  provided  by  CHL,  is  used  to  read  the  RMAIO  binary  file  and 
generate  a  3-D  ASCII  element  connection  file  named  R4ICR10.out.  This 
file  describes  the  juxtaposition  of  all  elements  and  nodes  in  the  RMAIO 
grid. 

c.  The  linkage  program  MAPPER.C  reads  R4ICR10.out  and  generates  the 
ICM  MAP  and  GEO  files.  The  number  and  position  of  ICM  cells  or 
boxes  are  determined  by  MAPPER  automatically.  Presently,  a  parameter 
statement  within  the  MAPPER  code  defines  “number_of_layers,”  which 
must  be  set  to  the  desired  number  of  layers.  This  parameter  can  be  used 
to  limit  the  output  to  a  single  layer  for  debug  purposes.  When  this 
parameter  is  set  to  the  correct  number  of  RMAIO  layers,  the  complete 
MAP  and  GEO  files  are  automatically  generated.  MAPPER  input  (i.e., 
R4ICR10.out)  and  output  files  for  a  simple  3  by  2  by  2  grid  are  presented 
in  Appendix  D. 

d.  The  time-varying  hydrodynamic  data  are  also  extracted  from  the  same 
RMAIO  binary  output  file.  This  is  accomplished  using  the 
FEMCONVT.f  program.  This  program  has  one  interactive  input,  which 
is  the  RMAIO  binary  file  name.  All  grid  parameters  are  read  from  the 
RMAIO  binary  file.  The  output  of  this  program  is  an  ASCII  ICM  HYD 
file.  At  this  time,  FEMCONVT.f  has  several  limitations: 
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•  ASCn  output  should  be  written  in  binary  to  reduce  space  requirements 
for  applications  with  large  output  files. 

•  Base  elevation  in  feet  relative  to  an  RMAIO  datum  must  be  hardwired. 

•  FEMCONVT  only  works  for  quadrilateral  elements. 

•  Parameter  statements  must  be  checked  to  ensure  the  code  will 
accommodate  the  size  of  the  grid. 

•  Vertical  diffusivities  required  by  ICM  are  not  computed  nor  output. 

•  There  is  no  provision  for  temporal  or  spatial  averaging  of  hydrodynamic 
output. 

e.  Two  additional  parameters  must  be  defined  in  the  ICM  file 

WQM_COM.INC.  NFEMHFFP  and  NFEMVFFP  are  the  maximum 
number  of  horizontal  and  vertical  flow  faces,  respectively. 
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3  Testing 


Volume  Balance  Testing 

Initial  testing  of  the  linkage  methodology  using  hydrodynamic  flow 
information  was  performed  on  a  simple  26  by  4  by  3  rectangular  contracted 
channel  comprised  of  quadrilaterals  as  shown  in  Figure  2.  To  test  both  local  and 
global  volume  conservation,  a  single  steady-state  flow  field  was  used.  Given  that 
the  flow  field  has  in  fact  converged  to  a  steady-state  solution,  successive 
applications  of  the  flow  field  using  the  steady-state  distribution  of  element 
volumes  as  an  initial  state  should  result  in  no  temporal  change  in  volume.  Using 
these  data,  ICM  was  run  using  a  1-hr  hydrodynamic  update  for  0.2  days,  or  two 
hydrodynamic  updates.  As  expected,  global  volume  conservation  was  achieved. 
However,  significant  variations  in  local  volume  conservation  were  observed 
throughout  the  grid.  This  problem  was  further  investigated  by  performing  flow 
balances  on  individual  cells.  In  regions  where  elements  were  orthogonal,  flow 
imbalances  on  the  order  of  1  percent  were  found,  where  the  error  in  flow  balance 
is  based  on  the  difference  in  the  sum  of  the  flows  scaled  by  the  maximum  flow  in 
or  out  of  that  cell.  In  regions  where  elements  were  nonorthogonal,  flow 
imbalances  were  on  the  order  of  60  percent.  Discussions  with  CHL  suggested 
that  this  problem  can  be  overcome  or  at  least  reduced  to  errors  on  the  order  of  a 
few  percent.  Means  of  correcting  the  flow/volume  imbalances  are  being 
pursued. 

Given  that  RMAIO  flow  fields  can  be  produced  with  small  but  possibly 
acceptable  volume  errors,  two  alternatives  can  be  adopted  that  will  result  in  a 
local  volume  balance.  The  simplest  option  is  to  use  the  flows  from  the  RMAIO 
hydrodynamic  processor,  but  not  use  the  resulting  volumes.  Under  this  scenario, 
time-varying  volumes  are  computed  within  ICM  using  the  previously  computed 
volume  and  the  input  flows,  thereby  ensuring  volume  and  mass  conservation. 

The  downside  of  this  option  is  that  if  a  persistent  net  loss  or  gain  is  realized  in 
one  or  more  cells,  problems  could  develop  during  long-term  simulations. 
Specifically,  cells  could  either  dry  out  or  become  unrealistically  large. 

Alternatively,  both  RMAIO  flows  and  volumes  can  be  used  by  ICM.  The 
RMAIO  flow  data  would  again  be  used  within  ICM  to  calculate  conserving 
volumes  for  each  hydrodynamic  update  interval.  The  difference  in  the  RMAIO 
volumes  and  ICM  volumes  can  then  be  used  to  determine  the  volume  errors  for 
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Figure  2.  Volume  balance  test  grid 

each  cell.  If  the  volume  error  is  ignored,  i.e.,  RMAIO  volumes  are  used,  the 
model  will  not  conserve  mass,  but  concentrations  will  be  conserved.  For 
example,  if  a  water  quality  constituent  of  10  ppm  is  introduced  into  a  bay  that 
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contains  10  ppm,  the  concentrations  will  all  remain  10  ppm,  but  some  mass  of 
the  constituent  will  be  gained  or  lost. 

Volume  and  mass  conservation  can  be  forced  when  using  both  RMAIO  flows 
and  volumes  by  introducing  a  volume  error  correction.  However,  under  this 
option  the  concentrations  will  not  remain  conservative,  or  10  ppm  everywhere  as 
in  the  previous  example.  The  error  correction  for  each  hydrodynamic  update 
interval  is  computed  and  added  (or  subtracted)  incrementally  over  the  ICM  time- 
steps  within  that  update  interval.  Thus,  the  ICM  volume  computed  at  the  end  of 
each  hydrodynamic  update  interval  will  be  identical  to  the  corresponding 
RMAIO  volume.  The  simplest  way  to  implement  this  correction  is  to  generate  a 
source  flow  by  dividing  the  total  volume  error  for  each  hydrodynamic  update 
interval  by  the  hydrodynamic  update  time  interval.  As  the  ICM  simulation 
proceeds,  a  volume  source/sink  can  be  added/subtracted  by  multiplying  the 
source  flow  by  the  ICM  time-step. 


Transport  Testing 

To  test  the  revisions  to  ICM  for  the  unstructured  grid  QUICKEST  multipliers, 
transport  testing  was  done  using  a  one-dimensional  rectangular  grid  set  up 
without  and  with  the  unstructured  grid  multipliers  as  derived  from  the  GEO  and 
MAP  files  generated  by  MAPPER.  The  RMAIO  model  was  not  used  for  these 
tests.  This  grid  consisted  of  25  boxes  100  m  in  length  with  1-m^cell  face  areas. 

A  constant  flow  rate  of  0.1  mVsec  was  specified  with  no  horizontal  diffusion. 
Specifying  a  constant  source  concentration  at  the  inflow,  the  predicted  square 
wave  concentrations  in  both  simulations  (i.e.,  with  structured  and  with 
unstructured  grid)  were  identical  as  shown  in  Figures  3  and  4.  The  algorithm 
was  tested  for  variable  grid  spacing  by  alternating  between  100-  and  50-m  box 
lengths  and  repeating  the  simulations.  Additionally,  both  positive  and  negative 
flows  were  tested.  All  runs  yielded  identical  results,  which  confirmed  correct 
implementation  of  QUICKEST. 

Additional  testing  of  transport  was  performed  using  a  steady-state  RMAIO 
hydrodynamic  flow  field.  A  rectangular  channel  RMAIO  grid  was  used 
consisting  of  10  elements  along  the  channel,  4  elements  across  the  channel,  and 
3  elements  deep  (see  Figure  5).  Output  from  RMAIO  was  used  to  drive  ICM. 
Exact  global  mass  conservation  was  verified  via  a  simulation  of  a  conservative 
tracer  spot  dump  and  simulation  of  a  uniform  tracer  concentration  throughout  the 
grid  and  boundaries.  In  addition,  a  constant  source  square  wave  test  was 
performed.  Despite  the  coarseness  of  the  grid,  the  expected  characteristics  of 
QUICKEST  are  preserved  in  that  minimal  undershoot  and  overshoot  and 
smearing  of  the  wave  occur  over  about  five  grid  points  (see  Figure  6). 
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QUICKEST  SQ.  WAVE  TEST 


Figure  3.  Square  wave  test  result  with  structured  grid  QUICKEST  multipliers 
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Figure  4.  Square  wave  test  result  with  unstructured  grid  QUICKEST  multipliers 
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4  Summary  and 
Recommendations 


A  linkage  methodology  has  been  developed,  implemented,  and  tested  that 
allows  CE-QUAL-ICM  to  use  unstructured  grid  hydrodynamic  information  from 
the  RMAIO  finite  element  model.  This  has  been  accomplished  through  the 
development  of  two  postprocessor  programs,  which  provide  flow  and  linkage 
information  to  ICM.  Initial  testing  of  volume  conservation  suggests  that  local 
volume  conservation  problems  inherent  within  RMAIO  limit  the  general  skill  of 
the  linkage  system;  however,  a  procedure  for  correcting  the  volume  errors  within 
ICM  is  outlined.  Irrespective,  results  of  the  present  study  suggest  that  further 
work  towards  eliminating  volume  conservation  errors  within  RMAIO  is 
warranted.  In  addition,  a  centroid-based  QUICKEST  advection  scheme 
consistent  with  the  unstructured  form  of  RMAIO  has  been  developed, 
implemented,  and  tested.  Testing  of  the  scheme  shows  that  it  is  mass 
conservative  and  that  it  provides  results  consistent  with  the  well-tested  structured 
grid  version  of  QUICKEST  previously  employed  in  ICM. 

With  respect  to  future  refinements  of  the  linkage  methodology,  in  addition  to 
either  correcting  for  or  eliminating  the  local  volume  conservation  errors  obtained 
from  RMAIO,  the  hydrodynamic  processor  FEMCONVT  needs  to  be  generalized 
to  accommodate  combinations  of  triangular  and  quadrilateral  elements. 
Furthermore,  additional  work  is  required  to  integrate  vertical  eddy  diffusivities 
within  RMAIO  SO  that  cell  face  values  can  be  provided  to  ICM.  Subsequent  to 
these  and  other  more  minor  improvements  outlined  in  Chapter  2,  the  complete 
linkage  methodology  should  be  tested  using  a  real  time-varying  field  application. 
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Appendix  A 

MAPPER  Description  and 
Source  Code  Listing 


The  following  is  a  brief  description  of  the  internal  steps  followed  by  the 
program  MAPPER.C  that  generates  the  MAP  and  GEO  files  needed  by  the  ICM 
model.  This  program  is  run  after  the  RMAIO  program  is  executed.  There  is  a 
program  called  R4ICR10.F  written  by  the  U.S.  Army  Engineer  Waterways 
Experiment  Station  Coastal  and  Hydraulics  Laboratory  that  extracts  the  two  files 
used  as  input  by  MAPPER.C  from  the  RMAIO  binary  output  file.  A  listing  of 
the  MAPPER.C  source  code  follows  the  description. 

The  step  by  step  operations  of  MAPPER.C  are  listed  below. 


Read  the  node-to-element  connection  table.  Allocates  each  “Element”  as  a 
column  and  assigns  an  ID. 

Subroutine  call:  readElementConnectionFile() 

Subroutine  description:  The  routine  opens  and  reads  the  element  to  node 
connection  table  ffle.  This  file  provides  information  that  uniquely  relates  an 
element  or  column  identification  number  to  the  associated  nodes.  Additional 
information  obtained  from  this  file  is  the  coordinates  of  each  node.  At  this  point, 
it  is  assumed  that  RMAIO  node  identification  numbers  are  being  used.  This 
routine  creates  a  unique  “arc”  for  each  side  of  the  surface  layer  elements,  which 
can  be  either  trapezoidal  or  triangular. 

Source  of  input  file:  Input  file  is  generated  by  R4ICR10.f. 


Read  boundary  information. 

Subroutine  call:  readBoundaryArcFile() 
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Subroutine  description:  Routine  reads  a  file  that  contains  information  about  all 
“special”  nodes  in  the  RMAIO  grid.  These  special  nodes  are  nodes  at  which 
some  physical  parameters  are  specified  such  as  head  values,  wall  boundaries, 
open  boundaries,  interior  boundaries,  etc.  This  information,  in  addition  to  the 
node  data  from  Step  1,  completely  describes  the  RMAIO  grid  structure. 

Source  of  input  file:  Input  file  is  generated  by  R4ICR10.f. 


Relate  each  arc  to  the  adjacent  columns. 

Subroutine  call:  mapColumnsToArcs() 

Subroutine  description:  Routine  determines  the  adjacent  column  identification 
numbers  for  each  arc. 


Relate  each  column  to  its  neighbors. 

Subroutine  call:  mapNeighborColumnsForEachColumnO 

Subroutine  description:  Routine  determines  for  each  column,  all  the  columns 
that  share  a  common  arc.  Thus,  a  list  is  generated  for  each  column  that  contains 
all  of  its  neighbor  columns. 


Determine  arc  type. 

Subroutine  call:  determineArcTypeFO 

Subroutine  description:  Routine  examines  each  arc  and  determines  its  type, 
which  is  defined  to  be  (a)  INTERNAL  -  arc  that  is  internal  to  the  grid, 

(b)  BOUNDARY  -  arc  lies  at  an  open  boundary,  and  (c)  WALL  -  arc  lies  at  the 
wall  of  the  grid.  The  arc  type  is  defined  using  the  nodal  boundary  data. 


Determine  upstream  and  downstream  columns. 

Subroutine  call:  mapUpstreamAndDownstreamColumnsForEachColumnF() 
Subroutine  description:  Routine  defines  adjacent  columns. 


Create  a  list  of  all  elements  in  the  3-D  grid. 

Subroutine  call:  createElementList() 

Subroutine  description:  Routine  will  generate  a  list  of  all  elements  in  the 
RMAIO  grid.  First,  it  generates  surface  layer  or  column  elements.  Next,  for 
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each  column,  it  will  define  elements  within  that  column  from  the  surface  down. 
The  order  of  elements  is  that  defined  by  RMAIO. 


Map  neighboring  elements. 

Subroutine  call:  mapElementToNeighborElements() 

Subroutine  description:  Routine  determines  neighboring  elements.  This  is 
accomplished  by  using  the  column  information  previously  defined. 


Check  for  internal  and  boundary  flow  faces. 

Subroutine  call:  generateHorizontalFlowFaceListF() 

Subroutine  description:  Routine  examines  each  arc,  its  type.  If  the  type  is 
“INTERNAL”  or  “BOUNDARY,”  it  wUI  create  a  flow-face  for  that  arc;  this  is 
done  in  layer  order  starting  a.t  the  surface  then  going  down.  For  each  defined 
flow  face,  the  adjacent  elements  of  that  flow  face  are  listed. 


Determine  far  upstream  elements  for  each  horizontal  flow  face. 

Subroutine  call:  mapHorizontaMowFaces() 

Subroutine  description:  Routine  determines  for  each  flow  face  its  far  upstream 
and  far  downstrejun  element  IDs  using  the  adjacent  element  ID  information. 


Create  the  vertical  flow  face  list. 

Subroutine  call:  mapYerticalFlowFacesQ 

Subroutine  description:  Routine  uses  the  coliunn  information  defined  above  and 
creates  a  vertical  flow  face  between  each  element.  It  also  does  the  vertical 
mapping  of  upstream  and  downstreeun  elements. 


Define  tags  that  relate  RMAIO  and  ICM  flow  faces. 

Subroutine  call:  setFlowFaceIdentifierTags() 

Subroutine  description:  Routine  examines  each  flow  face  and  tags  it  with  a 
unique  identifier.  This  unique  identifier  is  defined  as  (a)  horizontally  oriented 
flow  face  (the  top  center  node  ED  of  each  flow  face);  and  (b)  vertically  oriented 
flow  face  (the  RMAIO  element  ID  of  the  element  below  the  flow  face). 
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Number  ICM  flow  faces. 


Subroutine  call:  renumberAllFlowFacesForOutputF() 

Subroutine  description:  Routine  renumbers  the  flow  face  identifications 
consistent  with  the  ICM  convention. 

Calculate  surface  areas  and  the  centroid  locations. 

Subroutine  call:  calculateColumnSurfaceAreaAndCentroid() 

Subroutine  description:  Routine  will  calculate  surface  area  and  centroid  location 
using  the  surface  node  coordinates. 

Calculate  required  distances  between  centroids. 

Subroutine  call:  calculateDistanceBetweenCentroidsF() 

Subroutine  description:  Routine  calculates  the  following  length  scale 
information  for  each  flow  face: 

(a)  BIL  =  Distance  between  far  upstream  element  centroid  and  upstream 

element  centroid. 

(b)  BI  =  Distance  between  upstream  element  centroid  and  the 

downstream  element  centroid. 

(c)  BIR  =  Distance  between  downstream  element  centroid  and  far 

downstream  element  centroid. 

(d)  BID  =  Distance  between  flow  face  and  the  upstream  element  centroid. 

Print  “MAP  FILE”. 

Subroutine  call:  printFlowFaceListFQ 

Subroutine  description:  Routine  writes  “MAP”  file  data. 

Print  “GEO  FILE”. 

Subroutine  call:  printGeoFileFO 

Subroutine  description:  Routine  writes  “GEO”  file  data. 
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/*  PROGRAM  NAME:  MAPPER.c 
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#define  NUMBER_OF_LAYERS 


!#define  INTERNAL 
#define  BOUNDARY 
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int  bottomFlowFaceld; 

struct  ELEMENT_BLOCK_STRUCT  *elenientAbovePtr; 
struct  ELEMENT_BLOCK_STRUCT  *elementBelowPtr; 
}  ELEMENT_BLOCK; 


typedef  ELEMENT_BLOCK  *PELEMENT_BLOCK; 
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PCOLUMN_BLOCK  *lookupColumnPtrArray  =  NULL 


a< 

o 

o 

PQ 

a> 

T3 

O 

Z 

w 

Vi 

CQ 


H 
U 
D 
Di  , 

w  U 

o  ^ 

m 


w 

Q 

O 

:z: 


c/3 

<u 

T3 

U 

a< 


U 

O 

-I 

Q 

o 

z 

'W 

o 

3 


CU 

o 

5 

o 

•t3 

O 

•«-» 

X 

<0 

c 

4f- 

U 

D 

a:: 

(Z! 

o 

O 

h4 

03 

g'-b- 

is 

Z  o 

«  C 
o 

5  - 

*5  e 


m 

¥ 
■4^ 

CS 

c 

kl 

o 
o 
o 


U 

O 

J 

03 

03' 

Q 

O 


t3di 

u 

o 

03 

I 

w 

Q 

O 

12: 

CU 

■*• 

U 

o 

03 

03' 

Q 

O 

Z 

Cm 

U 

•o 

<u 

a. 

>4 


J 

J 

D 

Z. 

II 


o 

_o 

S 

(t> 

TJ 

O 

Z 

•o 

{3 

0) 

K 

a, 


U 

O 

J 

03 

03' 

Q 

O 

z 

CU 


c 

3 

o 

U 

(U 

•o 

o 

Z 

"S 

x> 

’w 


J 

J 

D 

Z 

II 

>» 

c3 


CL, 

O) 

■o 

o 

Z 

CL 

3 

O 

O 


U 

o 

J 

CQ 

I 

m 

Q 

o 

z 

CU 


CU  CU 

a  o 
_o  _o 

5  S 

<  < 


H 

U 

D 

DC 

H 

cc 


3 

J 

CL 

* 

H 

U 

D 


kSiC  oC 

o 

m 

lU 

u  " 

DC 
< 

o 
3 

Ui 
c/3 
Cm 
<U 

•o 
<u 
CL 
>4 


O 

m1 

03 

I 

U 

DC 

< 

o 

s 


X 

u 

CL 

* 

H 

U 

D 

DC 

H 

cn 

liC' 

U 

o 

m3 

«1 

< 


cn 


V 

TJ 

O 

a 


U 


o 

2  ^ 
Vi 


m3 

J 

D 

Z 


A8 


Appendix  A  MAPPER  Description  and  Source  Code  Listing 


typedef  ARC_BL0CK  *PARC_BL0CK; 


typedef  struct  FLOWFACE_BLOCK_STRUCT 
{  struct  FLOWFACE_BLOCK_STRUCT  *pLastFlowFaceBlockPtr; 
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PFLOWFACE_BLOCK  pHeadFlowFaceListBlockPtr  =  NULL; 

PFLOWFACE_BLOCK  pHeadSurfaceFlowFaceListBlockPtr  =  NULL; 

PH^OWFACE_BLOCK  pTailSurfaceFlowFaceListBlockPtr  =  NULL; 


PFL0WFACE_BL0CK  *lookupFlowFacePtrArray  =  NULL; 


/*  function  prototypes 


node_block 

*allocateNewNodeBlockF(); 
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PCOLUMN.BLOCK 
lookupColumnPtrF(int) ; 
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/*  MAIN  ROUTINE 


readElementConnectionFileO ; 
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mapElementToNeighborElementsO; 


printGeoFileFO; 
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if(pHeadNodeBlockPtr  ==  NULL) 

{  pHeadNodeBlockPtr  =  aNodeBlock; 


else 

{ lastNodePtr->nextNodeBlockPtr  =  aNodeBlock; 
aNodeBlock->lastNodeBlockPtr  =  lastNodePtr; 


rn  p3  Pfl 

c  c  c 

03  CO  C8 

3d  " 

§  H 

ass 

•T3  "O  D 

ea  ea 

£S1 

cii  cii  ca 


Appendix  A  MAPPER  Description  and  Source  Code  Listing 


A19 


retum(anElementBlockPtr); 


if(pHeadColumnListBlockPtr  ==  NULL) 

{  pHeadColumnListBlockPtr  =  aColunuiBlockPtr; 
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aFlowFaceBlockPtr->pLastFlowFaceBlockPtr  =  NULL; 
aFlowFaceBlockPtr->pNextFlowFaceBlockPtr  =  NULL; 
aFlowFaceBlockPtr->arcPtr  =  NULL; 


A22 


Appendix  A  MAPPER  Description  and  Source  Code  Listing 


void  calculateCentroidOfGeneralizedTriangle(x,y,cList) 
float  x[3],  y[3]; 
float  *cList; 


PCOLUMN_BLOCK  columnPtr; 
PARC_BLOCK  arcPtr; 
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/*  Each  column  has  4  potential  arcs  */ 


if(columnPtr->arcIds[i]  <  1)  continue; 
arcid  =  columnPtr->arcIds[i]; 
arcPtr  =  findArcPtrF(arcId); 
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for(i=l;  i<=nodeCount;  i++) 

{  nodePtr  =  findNodePtrF(nodeIds[i-l]); 


x[i]  =  nodePtr->coordinates[0]; 
y[i]  =  nodePtr->coordinates[l]; 
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columnPtr->centroidX  =  centroidX; 
columnPtr->centroidY  =  centroidY; 


columnPtr  =  columnPtr->pNextColumnBlockPtr; 
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upstreamColumnPtr  =  lookupColumnPtrF(upstreamColumnld); 


else 

upstreamColumnPtr  =  NULL; 
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retum(d); 


upstreamElementId  =  flowFacePtr->upstreamElementId; 
upstreamElementPtr  =  lookupElementPtrF(upstreamElementld); 
if(upstreamElementPtr  !=  NULL) 
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nodelPtr  =  findNodePtrF(nodel); 
xAl  =  nodelPtr->coordmates[0]; 
yAl  =  nodelPtr->coordinates[l]; 


node2Ptr  =  findNodePtrF(node2); 
xA2  =  node2Ptr->coordinates[0]; 
yA2  =  node2Ptr->coordinates[l]; 
ixl  =  (xAl  +  xA2)/2.0; 
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else 

lineType  =  INBETWEEN; 
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arcPtr  =  flowFacePtr->arcPtr; 
nodel  =  arcPtr->nodeId[0]; 
node2  =  arcPtr->nodeId[2]; 


nodelPtr  =  findNodePtrF(nodel); 
xAl  =  nodelPtr->coordinates[0]; 
yAl  =  nodelPtr->coordinates[l]; 
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yAl  =  nodelPtr->coordinates[l]; 


PELEMENT_BLOCK  upstreamElementPtr; 
int  upstreamElementld; 


cs 

X 


CO 

o 

5:3 


*0 


o 

^  3 

5  O 

53  T3 


U 

CO 


Appendix  A  MAPPER  Description  and  Source  Code  Listing 


A35 


farUpstreamElementld  =  flowFacePtr->farUpstreamElementId; 
farUpstreamElementPtr  =  lookupElementPtrF(farUpstreamElenientld); 
if(  (upstreamColumnPtr  !=  NULL)  &&  (farUpstreamElementPtr  !=  NL 


farUpstreamColumnld  =  farUpstreamElementPtr->columnId; 
farUpstreamColumnPtr  =  lookupColumnPtrFffarUpstreamCo 
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PELEMENT_BLOCK  farDownstreamElementPtr; 
int  farDownstreamElementld; 
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xl  =  downstreamColumnPtr->centroidX; 
yl  =  downstreamColunrinPtr->centroidY; 


x2  =  farDownstreamColumnPtr->centroidX; 
y2  =  farDownstreamColumnPtr->centroidY; 
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flowFacePtr  =  pHeadFlowFaceListBlockPtr; 


if(  (flowFacePtr->type  ==  INTERNAL)  II  (flowFacePtr->type  ==  BOUNDARY) ) 
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calculate_BID_F(flowFacePtr) ; 


flowFacePtr->distanceBetweenUpstreamCentroids  = 
calculate_BIL_F(flowFacePtr); 
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int  tempint; 


for(i=0;  i<3;  i++) 
for(j=i+l;j<3;j++) 
if(listl[i]>listl|j]) 

{  tempint  =  listl[i]; 
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PARC_BLOCK  checkIfArcExistsF(  int  list[3] ) 


NODE_BLOCK  *cuiTentNodePtr; 
int  i,  lineType; 
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if(  !strcmp(dummyStr,"GE") ) 
lineType  =  GE; 
else 


{  if(  !strcmp(dummyStr,"GNN") ) 
lineType  =  GNN; 
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currentColumnPtr->arcCount 
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cuiTentColumnPtr->arcIds[arcCount]  =  arcPtr->arcId 


index  =  0; 

for(arcCount=0;  arcCount<3;  arcCount++) 


PCOLUMN_BLOCK  columnPtr; 
PARC_BLOCK  arcPtr; 
int  arcid; 

int]  i,  iTemp; 
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else 

{ if(coluninPtr->columnId  >  arcPtr->adjacentColumn[l]) 

{  arcPtr->adjacentColunin[0]  =  arcPtr->adjacentColumn[l]; 
arcPtr->adjacentColumn[l]  =  columnPtr->colutnnId; 
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while(  boundaryNodePtr  !=  NULL) 

I  {  nodeld  =  arcPtr->nodeId[i]; 
if(boundaryNodePtr->boundaryNodeId  ==  nodeld) 
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while(arcPtr  !=  NULL) 

{  if((arcPtr->adjacentColumn[0]  !=  -1)  &&  (arcPtr->adjacentColumn[l]  !=-!)) 
arcPtr->type  =  INTERNAL; 


otherColumnId  =  searchColumnListForNamedArc(  columnPtr,  arcid  ); 
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retum(YES); 
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arcPtr2  =  findArcPtrF(arcId2); 


foundOppositeArc  = 
break; 
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arcPtr3  =  findArcPtrF(coluninPtr->arcIds[2]); 


{  arcPtr2  =  findArcPtrF(columnPtr->arcIds[0]); 
arcPtrS  =  findArcPtrF(columnPtr->arcIds[l]); 


U  I — II — I 

•o  o  — 


'  O  f— I 

0  0—1 


N<-4 

X) 

(D 

T3 

HH 

T3 

c 

Id 

CO 

CQ 

c 

■*-» 

W5 

o 

•w 

c« 

O 

O 

ij  _c 

C 

•o 

o 

C 

A 

S 

A 

Dh  ^ 
U  b; 

Vh 

T3 

Uh 

c 

A 

£ 

1 

<S 

i3 

CiH 

o 

cn 

U( 

£ 

o 

•o  o 
o  o 

i  X 

O 

O 

u 

A 

1 

O 

z 

73 

C 

o 

o 

o 

A 

f 

o 

o 

o 

A 

1 

o 

c3 

c3 

C  i 
«  CL, 

u 

IS 

i3 

U4 

11 

11 

£ 

II 

CL, 

cu 

II  « 

-  2 
Oi  C 

u 

0) 

li 

cs 

"O 

m 

T3 

l-H 

73 

O 

G 

CM 

•w 

cu 

T3 

O 

C 

T3 

O 

c 

o 

•a 

1> 

73 

^  II 

II 

o 

T3 

II 

It 

o 

o 

O 

o 

O 

cs 

<s 

C  X 


•o  13  OJ 
T3  'a  T3 


■O  T3  cn 
O  no  O 


Appendix  A  MAPPER  Description  and  Source  Code  Listing 


return(arcPtr2); 

else 

return(arcPtr3); 


return(NULL); 
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columnPtr->commonNeighbor2ColumnId[i]  = 


for(iCurrentCellLayerId  =  2; 
iCurrentCellLayerId<=iMaximumNumberOfLayers; 
I  iCurrentCellLayerId++) 
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elementPtr->elementId  =  iNumberOfCells; 


columnPtr  =  columnPtr->pNextColumnBlock;Ptr; 


c 

o 


0) 


=3 

fe  § 

■g  s 

C  D 

c  g 

o  S 

z  S 


-J  iS 

»  S 
'll  0^ 
^  CQ 

^  e 
^  l> 

^  B 

<L>  <l> 

m  « 

C  A 


"O 

•w 

c 

a> 

E 

(O 

U 

o 


nJ  c«  jr 

d  s  £ 

5  Gi-  w 

^  M  ^ 

il  -2  ^ 

*  +  "s 
%  Ji  I 
o  IE  V 


w 

X) 

< 


E  2 


1 

i3 

fi 

•£3 

Gm 

E 

U 

o 

Gm 

•w 

c 

13 

i 

o 

13 

2 

A 

A 

_3 

E 

■© 

o 

£ 

13 

a! 

4~> 

A 

II 

i! 

C 

1 

II 

w 

c 

1 

IX 

4-4 

c 

e 

<u 

o 

(D 

•w' 

GU 

c 

'S 

s-x 

G 

£ 

Ic 

£ 

0) 

^m4 

x: 

13 

(U 

& 

13 

II 


£ 


3 

O 

o 


Appendix  A  MAPPER  Description  and  Source  Code  Listing 


A61 


void 

mapElementToNeighborElementsO 


PCOLUMN_BLOCK  columnPtr,  otherColumnPtr,  farColumnPtr; 
irit  columnid,  otherColumnld,  farColumnId; 


c 


c 

E 

53 


a> 


« _ 

s 

3 

C 

B 

3 

6 

>< 

£  .2 
0)  — ’ 
c  c 


ii  o 

u  u 

Q  c 


u 

3 


d  £ 


u 

X) 

•y  S 

u  Z 
3  1 

I 

C 

3  X 
o  iS 


3  C 


c 

§ 

OJ 


lili 

U 

O 

pa 

h' 

z 

u 

s 

m 

j 

I? 

CL, 


O 

II 


3 


d  O 


-1 

J 

D 

Z 


1 

3 


S 

S 

‘S 


a 

c 

(U 

i: 

o 

o 


cd 


O 

1 

p 


X 

cd 


I! 

V 

•o 

HH 

U4 

o 

G 

g 

3 

o 


x: 

s: 


A62 


Appendix  A  MAPPER  Description  and  Source  Code  Listing 


columnPtr  =  pHeadColumnListBlockPtr; 


elementPtr  =  columnPtr->surfaceElementPtr; 
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flowFacePtr->adjacentColumn[0]  =  arcPtr->adjacentColumn[0] 
flowFacePtr->adjacentColumn[l]  =  arcPtr->adjacentColumn[l] 
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case  INTERNAL: 

flowFacePtr->type  =  INTERNAL; 
break; 
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PELEMENT.BLOCK  elementPtr; 
PFLOWFACE_BLOCK  flowFacePtr; 
PARC_BL0CK  arcPtr,  oppositeArcPtr; 


case  BOUNDARY: 

leftElementld  =  flowFacePtr->adjacentElement[0]; 
rightElementId  =  flowFacePtr->adjacentElement[l]; 
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flowFacePtr->farUpstreaniElementId  =  farLeftElementId; 


if(rightElementId  ==  -1) 
flowFacePtr->farDownstreaiTiElementId 
else 
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break; 


case  WALL: 
break; 
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UP!  IRREGARDLESS  of  other 
conventions 


cblumnPtr  =  pHeadColumnListBlockPtr; 
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arcPtr  =  arcPtr->pNextArcBlockPtr; 
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PCOLUMN_BLOCK 
findColuninPtr(int  columnld) 


Appendix  A  MAPPER  Description  and  Source  Code  Listing 


A77 


lookupColumnPtrF(int  columnld) 


if(lookupElementPtrArray  ==  NULL) 


* 
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lookupFlowFacePtrF(int  flowFaceld) 


PFLOWFACE_BLOCK  flowFacePtr; 


A80 
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PFLOWFACE_BLOCK  flowFacePtr; 
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else 

if(flowFacePtr->orientation  ==  VERTICAL) 


elementid  =flowFacePtr->upstreamElementId; 
elementPtr  =  lookupElementPtrF(elementld); 
flowFacePtr->identifierTag  =  elementPtr->rmalOElementId; 
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flowFacePtr  =  flowFacePtr->pNextFlowFaceBlockPtr; 
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case  WALL: 
faceType  = 
break; 


switch(flowFacePtr->orientation) 
{  case  HORIZONTAL: 
faceOrientation  =  'H'; 
break; 
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PELEMENT.BLOCK  elementPtr; 
PFL0WFACE_BL0CK  flowFacePtr; 


PELEMENT_BLOCK  upstreamElementPtr; 
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retum(intList); 


void 

printFlowFaceListFO 
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while(elementPtr  !=  NULL) 
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modelFlowFaceld  =  flowFacePtr->modelFlowFaceId; 


A88 


< 


c 

u 


o 

II 

T3 

C 

i 

0) 


a 

Vi 

CU 


CS 

(U 

■*-* 

CO 

Oh 


U 

Oh 


ts 

Oh 

U 

o 

es 

Oh 

o 

C 


Oi 

< 

Q 

Z 

P 

o 

m 

II 


V 

Oh 


Oh 

<1> 

O 

r'* 

Oh 

O 


I 

VO 


o 


T3 

VO 


Appendix  A  MAPPER  Description  and  Source  Code  Listing 


faceType, 

faceOrientation, 

flowFacePtr->identifierTag); 
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coluninPtr  =  columnPtr->pNextColumnBlockPtr; 


/*  Write  out  the  vertical 
flow  faces  in  a  column  info 


* 
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columnPtr  =  columnPtr->pNextColuninBlockPtr; 
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columnPtr  =  pHeadColumnListBlockPtr; 


elementPtr  =  columnPtr->surfaceElementPtr; 
surfaceElementPtr  =  elementPtr; 
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fprintf(outf,"%6d  \n",elementAboveId); 


elementPtr  =  elementPtr->pNextElementBlockPtr; 
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BI  =  flowFacePtr->distanceBetweenCentroids; 

BID  =  flowFacePtr->downstreamDistanceBetweenCentroidAndFace; 
BIR  =  flowFacePtr->distanceBetweenDownstreamCentroids; 


infile  =  fopen("r4icrl0.output","r"); 
if(infile  ==  NULL) 

{  printfC'Unable  to  open  'r4icrl0.output'\n"): 
fflush(stdout); 
exit(O): 
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surfaceNodeList[surfaceNodeId].ndep  =  ndep; 
surfaceNodeList[surfaceNodeId].nref  =  nref; 


for(i=0;  i<3;  i++) 
fgets(lineBuffer,  132, infile); 
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TO  SOMETHING  MORE  GENERAL 


open(unit=  1 0,file='femconvT.out', status  ='UNKNOWN') 
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nnpe(j)=l 

endif 

CONTINUE 


LOOP  THROUGH  ALL  THE  2D  SURFACE  ELEMENTS  IN  THE  GEOMETRY 
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IDND(IDEL,4)  =  NREF(IDEL)+1 
IDND(IDEL,5)  =  NREF(N0DR)+2 
IDND(IDEL,6)  =  NREF(N0DR)+1 
IDND(IDEL,7)  =  NODR 
IDND(IDEL,8)  =  JNN 


NNPE(IDEL)  =  8 
ELSE 

IDEL=NREF(JNN)+J-1 

JNM=JNM+1 

LIS(JNM)=IDEL 
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IDND(IDEL,4)  =  NREF(IDEL)+1 
IDND(IDEL,5)  =  NODR 
IDND(IDEL,6)  =  IDEL 
NNPE(IDEL)  =  6 
ELSE 
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IDND(IDEL,6)  =  IDEL 
NNPE(IDEL)  =  6 
ENDIF 
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Now  calculate  fluxes 


THIS  IS  CURRENTLY  SETUP  FOR  QUADRILATERAL  ONLY 
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V4=XVEL(2,IDND(IFACE,4)) 
U5=XVEL(  1  ,IDND(IFACE,5)) 
V5=XVEL(2,IDND(IFACE,5)) 
H5=XVEL(3,IDND(IFACE,5)) 
U6=XVEL(1,IDND(IFACE,6)) 
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eg 


temporary  variables  for  jacobians 
X1Z1=X1  *Z1 
X3Z1  =  X3  *  Z1 


X5Z1  =  X5  *  Z1 
X1Z3  =  X1  *Z3 
X3Z3  =  X3  *  Z3 
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FY  =  (V1*(-X1Z3+X3Z1-2.*X3Z5+2.*X5Z3) 
+2.*V2*(-4.*X3Z1+X5Z1+4.*X1Z3 
-3.*X5Z3-X1Z5+3.*X3Z5) 
+V3*(X3Z1+X5Z1-X1Z3+X5Z3-X1Z5-X3Z5) 
+2.*V4*(-3.*X3Z1+X5Z1+3.*X1Z3 
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JNV  =  JNV+1 
LISE(JNV)=IE 
IB1=N0P(IE,1) 
IB3=NOP(IE,3) 


IB5=NOP(IE,5) 

IB7=NOP(IE,7) 

IT1=N0P(IE,13) 

IT2=NOP(IE,14) 

IT3=NOP(IE,15) 
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ZT7=CORD(IT7,3) 

ZB  1=(ZB  1-A0(IB  1))/(ELEV-A0(IB  D) 
*XVEL(3,IB1)+A0(IB1) 
ZB3=(ZB3-AO(IB3))/(ELEV-AO(IB3)) 


*XVEL(3,IB3)+AO(IB3) 

ZB5=(ZB5-AO(IB5))/(ELEV-AO(IB5)) 

*XVEL(3,IB5)+AO(IB5) 

ZB7=(ZB3-AO(IB7))/(ELEV-AO(IB7)) 

*XVEL(3,IB7)+AO(IB7) 
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+X7Y3-X3Y5+X7Y5+2.*XlY7-X3Y7-X5Y7)/i: 
+(ZB3-ZT3)*(2.*X3Y1-X5Y1-X7Y1-2.*X1Y3+ 
2.  *X5  Y3+X 1 Y5-2.  *X3  Y5+X7  Y5+X 1 Y7-X5  Y7)) 
+(ZB5-ZT5)*(X3Y1-X7Y1-X1Y3+2.*X5Y3-X7^ 
-2.  *X3  Y5+2.  *X7  Y5+X 1 Y7+X3  Y7-2.  *X5  Y7)/ 12 


+(ZB7-ZT7)*(X3Y1+X5Y1-2.*X7Y1-X1Y3+X5Y3 
-X 1 Y5-X3  Y5+2.  *X7  Y5+2.  *X  1 Y7-2.  *X5  Y7)/12. 
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(-X3Y1+X7Y1+X1Y3-2.*X5Y3+X7Y3+2.*X3Y5- 

2.*X7Y5-XlY7-X3Y7+2.*X5Y7)/36.) 

TMP6=W6*((2.*X3Y1+X5Y1-3.*X7Y1-2.*X1Y3+3.*X5Y3- 

X7Y3-X1Y5-3.*X3Y5+4.*X7Y5+3.*X1Y7+ 

X3Y7-4.*X5Y7)/36.+ 
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write( 10,5000)  JNV 
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Several  modifications  were  made  to  the  original  ICM 
source  code  to  allow  it  to  be  driven  by  the  RMAIO 
hydrodynamic  model.  The  changes  made  are  detailed  below: 
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HHtagToItagMap(NFEMHFFP),  HVtagToItagMap(NFEMVFFP), 
FemCellToIcmCell(NBP) 

common  /lookup/  HHtag,  HVtag,  Itag,  HHtagToItagMap, 
HVtagToItagMap,  FemCellToIcmCell 
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2020  format(/) 

2030  format(I8,9I8) 


open  (geo,file=GEOFN,status='OLD') 
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1030  format(el5.9) 
1040  format(/) 


*****  Time-invariant  hydrodynamic  data 
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logical  function  initializeHydro() 
include  'model.inc' 


integer  HHtag(NFEMHFFP),  HVtag(NFEMVFFP),  Itag(NQFP), 

HHtagToItagMap(NFEMHFFP),  HVtagToItagMap(NFEMVFFP), 
FemCellToIcmCell(NBP) 
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do  i=l,NFEMHFFP 
HHtagToItagMap(i)  = 


unusedVFaceCount 
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HVtagToItagMap(HVtagValue)  =  Itagindex 

if(ltaglndex  .eq.  0)  unusedVFaceCount  =  unusedVFaceCount+1 


Modification  4: 

A  subroutine  for  reading  the  time- variant  hydro  data  was  written 
to  read  the  RMAIO  hydro  output  data.  This  routine  replaces  ICM's 
HYDRO  subroutine.  New  subroutine  code  follows: 
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NOTE  the  PARAMETER  NFEMFFP  MEANS  .. 
Number  of  Finite  Element  Model  Flow 
Faces  Parameter 

NFEMFFP  =  nHtagsl  +  nHtags2 


icm_box_id:  id  of  box  in  ICM 
fem_box_id;  id  of  box  in  RMAO 
volume:  box  volume 
top_area:  area  of  box  top 
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unusedVFaceCount  =  0 
read(hyd,100,end=50)  nHtags2 
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*********  Initialize  volumes 
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Modification  7: 

The  mainline  initialization  of  the  BL  array  was  eliminated: 


*********  Initialize  box  lengths 
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LOGICAL  END_OF_FILE,  readHydro 
DIMENSION  MASS(0:NBP,NCP) 


10000  IF  (NWQMR.GE.NHMR)  THEN 
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IF  (CONSERVE_MASS)  THEN 


DO  10047  B=1,NB 
DO  10045  JC=1,NAC 
MASS(B,AC(JC))  =  C1(B,AC(JC))*V1(B) 
10045  CONTINUE 
V1(B)  =V2(B) 
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BENFLXDP 

BENFLXPC 
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IWCMP  =  CMASS(15)+CMASS(16)+CMASS(17)+CMASS(18) 
IWCMC  =  CMASS(4)+CMASS(5)+CMASS(6)+CMASS(7)+CMASS(8) 
+CMASS(9) 

IWCMS  =  ASCD*CMASS(5)+CMASS(21)+CMASS(22) 

DO  10058  BB=1,NBB 


IF  (.NOT.FLOW)  THEN 
DO  10100  F=1,NQF 
Q(F)  =  0.0 
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DO  10140  F=NHQF+1,NQF 
DIFF(F)  =  ZDFMUL*DIFF(F) 
10140  CONTINUE 


*******  ASCII  input  FORMAT  statements 
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ZD(i)  is  depth  below  surface 
of  top  of  box  i! 
cell  depth  below  equals 
cell  depth  above  +  thickness 
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RMA10_HYDRO 


*****  Horizontal  advection  and  diffusion  multipliers 
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IF  (LEFT_FLOWB(F))  BL(IB(F),QD(F))  =  BL(JB(F),QD(F)) 
IF  (RIGHT_FLOWB(F))  THEN 


BL(JB(F),QD(F))  =  BL(IB(F),QD(F)) 

BL(JRB(F),QD(F))  =  BL(IB(F),QD(F)) 

END  IF 

IF  (RIGHTP1_B0UNDARY(F))  BL(JRB(F),QD(F))  =  BL(JB(F),QD(F)) 
SFl(F)  =  MIN(BL(IB(F),QD(F)),BL(JB(F),QD(F))) 


BI(F)  =2.0*BID(F) 
BIR(F)  =  BI(F) 

END  IF 

SFl(F)  =BI(F) 
SF2(F,1)  =BI(F)**2 
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TP2(F,2)  =  (2.0*BID(F)-BI(F)-BIR(F))*SF1(F) 
TP3(F,2)  =  (2.0*BID(F)-BI(F))*SF1(F) 


Appendix  D 

MAPPER  Input  and  Output 
Files  for  Example  Grid 


This  appendix  presents  an  example  grid  and  the  associated  input  and  output 
files  for  MAPPER.C.  The  example  consists  of  a  3  by  2  by  2  grid  domain  of 
12.0  m^,  with  element  dimensions  1.0  by  1.0  by  1.0  hl  The  node  numbers  are 
shown  at  the  comers  and  facial  mid-points,  and  the  element  numbers  are 
indicated  in  the  center  of  each  element  in  Figure  Dl.  Plan  views  are  shown  in 
Figure  Dl  for  the  two  layers  at  the  top,  middle,  and  bottom  of  each  layer.  The 
node-element  connectivity  output  file  from  RC4IC10.f  is  shown  for  the  example 
grid  in  Table  Dl.  This  file  provides  the  input  required  for  MAPPER.C.  The 
output  files  from  MAPPER.C,  MAP  and  GEO,  are  shown  for  the  example  grid  in 
Tables  D2  and  D3,  respectively.  Following  is  a  description  of  the  steps  to 
conduct  this  example,  along  with  explanations  for  the  input  and  output  files. 


MAPPER  Input  File  Description 

The  hydrodynamic  grid  is  initially  generated  by  the  FastTABS  program, 
which  is  a  preprocessor  and  postprocessor  for  two-dimensional  (2-D)  finite 
element  models  maintained  by  the  U.S.  Army  Engineer  Waterways  Experiment 
Station  Coastal  and  Hydarulics  Laboratory.  The  user  interactively  constracts  a 
grid  within  FastTABS  that  then  outputs  a  mesh  geometry  file  for  later  input  into 
RMA.  All  information  contained  in  this  file  is  subsequently  written  by  RMA  to 
what  is  referred  to  as  the  RMA  binary  output  file.  This  RMA  output  file 
ultimately  wUl  contain  a  complete  geometric  description  of  the  grid  in  3-D  along 
with  the  time- varying  hydrodynamic  flow  fields.  The  RC4IC10.f  postprocessing 
program  reads  the  RMA  binary  output  file  and  creates  the  node-element 
coimectivity  file,  which  is  the  input  file  for  MAPPER.C 

The  node-element  connectivity  file  (i.e..  Table  Dl)  is  described  as  follows. 
Information  in  the  first  section  of  this  file  describes  the  2-D  mesh  geometry  of 
the  grid.  Each  line  begins  with  a  line  type  code  consisting  of  two  or  three  letters. 
Code  definitions  follow  below. 
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Table  D.l.  Node-Element  Connectivity  File  for  Example  Grid  from  RC4IC10.f  Output  File 
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Table  D.2.  MAP  File  Generated  for  Example  Grid 
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Table  D.3.  GEO  File  Generated  for  Example  Grid 
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T1  thru  T3:  Titles.  T1  and  T2  are  optional,  T3  is  necessary. 

SI:  Signifies  English  or  SI  units.  English  =  0,  SI  =  1. 

$L:  logical  unit  numbers  used  by  RMA  specific  routines. 

GO:  Used  to  input  the  nodes  of  a  GO  string  (RMA  term).  Used  as  a 
reference  line  for  a  renumbering  of  node  identifier  numbers  done  by  RMA 
internally. 

GE:  Element  definitions. 

Field  Description 

I  The  first  integer  is  the  element  identifier. 

2-9  The  next  eight  integers  define  the  nodal  coimectivity  of  the  element. 
The  node  ordering  is  coimter-clockwise  beginning  at  a  comer  then 
sequentially  around  the  perimeter.  Eight  integers  are  always  used. 
Triangular  elements  that  have  only  six  nodes  will  have  the  remaining  two 
integers  set  to  zero. 

10  Material  identifier  (internal  to  RMA) 

II  A  real  value.  Represents  the  direction  of  the  eddy  viscosity  tensor 
for  the  element  and  is  usually  left  as  zero. 

GNN:  Comer  node  coordinates. 

Field  Description 

1  Integer.  Node  identifier. 

2  Real.  X  coordinate. 

3  Real.  Y  coordinate. 

4  Real.  Z  coordinate. 

Midside  nodes  are  not  listed.  They  can  be  calculated  from  comer  nodes. 

Information  in  the  second  section  Table  D1  describes  the  depth  dimension  of 
the  grid.  This  section  begins  with  the  specification  of  the  following  parameters: 

NP  -  Total  number  of  nodes  in  grid. 

NE  -  Total  number  of  elements  in  grid. 

NPM  -  Total  number  of  nodes  in  the  surface  plane. 
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NEM  -  Total  number  of  elements  in  surface  plane. 

Next,  the  number  of  nodes  beneath  each  surface  node,  including  the  surface 
node,  and  an  integer  indicating  the  numbering  of  this  column  of  nodes  are 
specified. 

NDEP  -  Total  number  of  nodes  in  a  column.  A  column  exists  for  each  node 
in  the  surface  plane. 

NREF  -  Nodes  in  the  surface  plane  are  numbered  sequentially.  However,  all 
other  nodes  are  numbered  in  sequential  order  from  the  node  beneath  the 
surface  node  to  the  bottom  within  a  node  column.  This  integer  is  the  starting 
node  identifier  for  the  first  node  beneath  the  surface  plane  node  in  a  node 
column  plus  1. 

The  final  section  of  Table  D1  indicates  for  each  node  within  a  node  column 
whether  it  is  internal  to  the  grid  or  positioned  on  the  boundary.  The  following 
definitions  apply. 

NSURF:  An  integer  associated  with  each  node  column  indicating  the  type  of 
node  column.  Possible  values  are  as  follows: 

0  -  Internal  node  column. 

1000  -  Wall  boundary  exists  at  this  node  column. 

1200  -  Head  boundaiy  exists  at  this  node  column. 

31000  -  M2  tidal  boundary  exists  at  this  node  column. 


MAPPER  Output  Files  Descriptions 

The  program  MAPPER.  C  reads  the  node-element  connectivity  file  and 
generates  the  output  files  MAP  and  GEO  files  used  as  input  by  ICM.  When 
MAPPER  reads  the  node-element  connectivity  file,  it  builds  an  internal  RMA 
grid.  MAPPER  is  aware  of  how  RMA  numbers  its  nodes  and  elements.  After 
completing  the  internal  RMA  grid,  MAPPER  can  then  write  out  the  linkage  and 
geometry  data  ICM  needs.  MAPPER  will  create  the  MAP  and  GEO  files.  The 
MAP  file  contains  the  RMA  element  to  ICM  cell  correspondence  and  flow  face 
defimtions.  The  GEO  file  specifies  the  vertical  relationship  between  ICM  cells 
and  lists  the  QUICKEST  weighting  distances. 


ICM  MAP  File  Description 

The  MAP  file  contains  four  sections  of  data.  A  brief  explanation  of  each 
section  follows  below. 
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The  first  section  lists  the  correspondence  between  ICM  cell  identifiers  and 
RMA  element  identifiers.  Note  that  the  two  models  have  different  numbering 
schemes  for  their  cells.  This  section  establishes  the  relationship  between  the  two 
numbering  schemes.  The  FEMCONVT.F  program,  which  reads  the  RMA  binary 
output  file,  writes  out  data  in  terms  of  RMA  internal  element  numbers  and  node 
numbers;  therefore,  the  information  in  this  section  will  allow  a  unique 
identification  of  each  cell  and  flow  face. 

RMA  numbers  elements  sequentially  in  the  order  of  their  creation  for  the 
surface  layer.  Next,  it  numbers  its  elements  from  the  bottom  up  to  the  layer 
beneath  the  surface  layer  going  from  one  column  to  the  next  in  the  order  of  the 
element  numbering  in  the  surface  layer.  The  ICM  numbering  scheme 
sequentially  numbers  the  cells  within  a  layer  starting  with  the  surface  layer  then 
progressing  to  the  layer  beneath  till  the  last  layer  has  been  processed. 

Definitions  of  the  parameters  for  this  section  foUow. 

ICM  -  Unique  ICM  cell  identifier. 

FEM  -  Unique  RMA  cell  identifier  corresponding  to  the  appropriate  ICM 
cell. 

The  second  section  identifies  each  flow  face  in  the  grid.  Definitions  of  the 
parameters  for  flow  face  definition  follow. 

FFBD  -  Unique  flow  face  identifier. 

ILB  -  QUICKEST  far  upstream  weighting  cell  identifier. 

IB  -  Cell  identifier  for  flow  donor  cell. 

JB  -  Cell  identifier  for  flow  receiver  cell. 

JRB  -  QUICKEST  far  downstream  weighting  cell  identifier. 

TYPE  -  Character  variable  indicating  type  of  flow  face.  B  indicates 
Boundary.  I  indicates  Internal  flow  face. 

ORIENT  -  Character  variable  indicating  direction  orientation.  H  indicates 
face  is  oriented  for  horizontal  flow.  V  indicates  face  is  oriented  for  vertical 
flow. 

LAYER  -  Indicates  what  layer  the  flow  face  is  in. 

The  third  section  lists  the  number  of  vertical  flow  faces  in  a  cell  column 
beneath  each  surface  cell.  Definitions  of  the  parameters  for  this  section  follow. 

SB  -  Surface  cell  (box)  identifier. 

NVF  -  Number  of  vertical  flow  faces  beneath  the  surface  cell. 
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The  fourth  section  lists  the  vertical  flow  face  identifiers  for  each  ICM  cell 
column  iterated  in  the  third  section.  Definitions  of  the  parameters  necessary  for 
this  section  follow. 

SB  -  Surface  cell  (box)  identifier. 

VFN  -  List  of  all  flow  face  identifiers  oriented  vertically  beneath  a  surface 
cell. 


ICM  GEO  File  Description 

The  GEO  file  contains  tree  sections  of  data.  A  brief  explanation  of  each 
section  follows  below. 

The  first  section  identifies  for  each  surface  cell  the  bottom  cell  directly 
beneath  it.  Definitions  of  the  parameters  for  this  section  follow. 

SB  -  Surface  cell  (box)  identifier. 

BBN  -  Bottom  cell  identifier. 

The  second  section  identifies  for  each  cell,  in  increasing  order,  the  cell 
directly  above  it.  Definition  of  the  parameter  for  this  section  follows. 

BU  -  Identifier  of  cell  above.  Cell  identifier  is  implicit  in  line  order. 

The  third  section  specifies  the  QUICKEST  weighting  distances  for  each  flow 
face.  Definitions  of  the  parameters  for  this  section  follow. 

BI  -  Centroid  to  centroid  distance  between  IB  cell  JB  cell. 

BID  -  Distance  from  centroid  of  IB  cell  to  flow  face. 

BIL  -  Centroid  to  centroid  distance  between  ILB  cell  and  IB  cell. 

BIR  -  Centroid  to  centroid  distance  between  JB  cell  and  JRB  cell. 

TYPE  -  Character  variable  indicating  type  of  flow  face.  B  indicates 

Boundary.  I  indicates  Internal  flow  face. 

LAYER  -  Indicates  what  layer  the  flow  face  is  in. 
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for  reducing  this  burden,  to  Washington  Heacfouarters  Services,  Directorate  for  Irtfoimation  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the 
Office  of  Mar»gement  and  Buciget,  Paperwork  Reduction  Project  (0704-0186),  Washington,  DC  20503. 
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